function [S, ST] = ek_struct(A, cholesky)
%BUILD_RK_STRUCT Build a struct for building rational Krylov subspaces.
%

if issparse(A)
    if exist('cholesky', 'var') && cholesky
        [R, g, pA] = chol(A, 'vector');
        
        if (g ~= 0)
            % warning('Matrix is not posdef, resorting to LU');
            [S, ST] = ek_struct(A, false);
            return;
        end
        
        qA(pA) = 1:size(A,1);
        
        S = struct(...
            'solve', @(nu, mu, x) sparse_solve(nu, mu, R', R, pA, qA, x), ...
            'multiply', @(rho, eta, x) rho * A * x - eta * x, ...
            'isreal', isreal(A), ...
            'nrm', normest(A, 1e-2));
        ST = S;
        
    else
        [LA, UA, pA, qA] = lu(A, 'vector');
        
        iqA(qA) = 1:size(A,1);
        ipA(pA) = 1:size(A,1);        
        
        S = struct(...
            'solve', @(nu, mu, x) sparse_solve(nu, mu, LA, UA, pA, iqA, x), ...
            'multiply', @(rho, eta, x) rho * A * x - eta * x, ...
            'isreal', isreal(A), ...
            'nrm', normest(A, 1e-2));
        
        if nargout >= 2
            ST = struct(...
                'solve', @(nu, mu, x) sparse_solve(nu, mu, UA', LA', qA, ipA, x), ...
                'multiply', @(rho, eta, x) rho * A' * x - eta * x, ...
                'isreal', S.isreal, ...
                'nrm', S.nrm);
        end
    end
    
elseif isa(A, 'hodlr')
    [LA,UA] = lu(A);
    
    S = struct(...
        ... % 'solve', @(nu, mu, x) (nu * A - mu * eye(size(A), 'like', A)) \ x, ...
        'solve', @(nu, mu, x) lu_solve(nu, mu, LA, UA, x),...
        'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
        'isreal', isreal(A), ...
        'nrm', norm(A));
    
    if nargout >= 2
        ST = struct(...
            ...% 'solve', @(nu, mu, x) (nu * A' - mu * eye(size(A), 'like', A)) \ x, ...
            'solve', @(nu, mu, x) lu_solve(nu, mu, UA', LA', x), ...
            'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
            'isreal', S.isreal, ...
            'nrm', S.nrm);
    end
    
elseif isa(A, 'hss')
    ULV = ulv(A);
    
    S = struct(...
        ... % 'solve', @(nu, mu, x) (nu * A - mu * eye(size(A), 'like', A)) \ x, ...
        'solve', @(nu, mu, x) ulv_solve_wrapper(nu, mu, ULV, x),...
        'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
        'isreal', isreal(A), ...
        'nrm', norm(A));
    
    if nargout >= 2
        ULVT = ulv(A');
        ST = struct(...
            ...% 'solve', @(nu, mu, x) (nu * A' - mu * eye(size(A), 'like', A)) \ x, ...
            'solve', @(nu, mu, x) ulv_solve_wrapper(nu, mu, ULVT, x), ...
            'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
            'isreal', S.isreal, ...
            'nrm', S.nrm);
    end
else
    % In this case we are probably dealing with a dense matrix, so we try
    % Cholesky first in the symmetric case, and otherwise resort to a
    % standard LU factorization
    if exist('cholesky', 'var') && cholesky
        [R, g] = chol(A);
        pA = (1 : length(A));
        
        if (g ~= 0)
            % warning('Matrix is not posdef, resorting to LU');
            [S, ST] = ek_struct(A, false);
            return;
        end
        
        S = struct(...
            'solve', @(nu, mu, x) sparse_solve(nu, mu, R', R, pA', pA, x), ...
            'multiply', @(rho, eta, x) rho * A * x - eta * x, ...
            'isreal', isreal(A), ...
            'nrm', normest(A, 1e-2));
        ST = S;
        
    else
        % We use this syntax to be able to call sparse_solve even if this
        % is the case of a generic matrix.
        [LA, UA, pA] = lu(A, 'vector');
        qA = 1 : length(A); 
        
        S = struct(...
            'solve', @(nu, mu, x) sparse_solve(nu, mu, LA, UA, pA, qA, x), ...
            'multiply', @(rho, eta, x) rho * A * x - eta * x, ...
            'isreal', isreal(A), ...
            'nrm', normest(A, 1e-2));
        
        if nargout >= 2
            ST = struct(...
                'solve', @(nu, mu, x) sparse_solve(nu, mu, UA', LA', pA', qA', x), ...
                'multiply', @(rho, eta, x) rho * A' * x - eta * x, ...
                'isreal', S.isreal, ...
                'nrm', S.nrm);
        end
    end    
end

end

